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A Module for Radiation Hydrodynamic Calculations With ZEUS-2D Using 

Flux-Limited Diffusion 

N. J. Turner^ & J. M. Stone^ 
ABSTRACT 

A module for the ZEUS-2D code is described which may be used to solve the equa- 
tions of radiation hydrodynamics to order unity in vjc^ in the flux-limited difl'usion 
(FLD) approximation. In this approximation, the tensor Eddington factor f which 
closes the radiation moment equations is chosen to be an empirical function of the 
radiation energy density. This is easier to implement and faster than full-transport 
techniques, in which f is computed by solving the transfer equation. However, FLD is 
less accurate when the flux has a component perpendicular to the gradient in radiation 
energy density, and in optically thin regions when the radiation fleld depends strongly 
on angle. 

The material component of the fluid is here assumed to be in local thermodynamic 
equilibrium. The energy equations are operator-split, with transport terms, radiation 
diffusion term, and other source terms evolved separately. Transport terms are applied 
using the same consistent transport algorithm as in ZEUS-2D. The radiation diffu- 
sion term is updated using an alternating-direction implicit method with convergence 
checking. Remaining source terms are advanced together implicitly using numerical 
root-finding. However when absorption opacity is zero, accuracy is improved by instead 
treating the compression and expansion source terms using a time-centered differencing 
scheme. 

Results are discussed for test problems including radiation-damped linear waves, 
radiation fronts propagating in optically-thin media, subcritical and supercritical radi- 
ating shocks, and an optically-thick shock in which radiation dominates downstream 
pressure. 

Subject headings: hydrodynamics — methods: numerical — radiative transfer 



1. INTRODUCTION 

In many astrophysical systems, radiation is the dominant energy transport mechanism. In some 
circumstances, including winds from Wolf-Rayet stars (Lucy & Abbott 1993) and asymptotic giant 
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branch stars (Habing 1996), interiors of massive stars (Kippenhahn & Weigert 1990), supernova 
blast waves (Arnett et al. 1989), and compact object accretion flows (Klein et al. 1996) and 
disks (Shakura & Sunyaev 1973), radiation comprises a substantial fraction of the total energy 
density, momentum density, or pressure. The dynamics of such flows must be modeled using a 
radiation hydrodynamical (RHD) approach, in which energy and momentum conservation laws for 
the radiation field are solved along with those for the gas. The outcomes of these calculations 
depend on the exchange of internal energy and momentum between material and radiation. The 
length and time scales associated with the exchanges are determined by atomic processes, and may 
be orders of magnitude shorter or faster than dynamical scales, so that implicit time differencing 
is usually necessary. An accurate description of the angular dependence of the radiation field may 
be important when parts of the gas are optically thin, localized sources of radiation are present, or 
shadows are cast within the flow. It is possible to compute the angular dependence by solving the 
transfer equation along rays through the points of interest. However, the number of rays needed 
to adequately specify the angular dependence at every point is often large. Such a method has 
been implemented in two spatial dimensions in the ZEUS-2D magnctohydrodynaniics (MHD) code 
(Stone & Norman 1992b) by Stone, Mihalas, & Norman (1992), but is complex and computationally 
intensive. Despite the promise of this full-transport method, it must be considered still under 
development. In this paper we describe a method for calculations using ZEUS-2D, in which the 
angular dependence of the radiation fleld is assumed to be given by the flux-limited diffusion (FLD) 
approximation. This method is simple, robust, and relatively cheap in computer time. We assume 
the gas is in local thermodynamic equilibrium (LTE) at a temperature which need not correspond 
to that of the radiation field. Frequency dependence of the opacities is neglected, though it may 
be included in a straightforward fashion. The RHD equations solved are discussed in §2, the FLD 
approximation in §3, and the numerical algorithm in §4. Results of test calculations are presented 
in §5, and the advantages and limitations of the method are summarized in §6. 



Much as the gas dynamical equations are derived by taking velocity moments of the material 
particle kinetic equation, or Boltzmann equation, the conservation laws for the radiation field 
are generated by taking angular moments of the photon kinetic equation, or radiation transport 
equation. In a frame comoving with the radiating fluid, assuming LTE, and to order unity uiv/c, 
the coupled equations of RHD are (Mihalas k, Mihalas 1984) 



2. EQUATIONS OF RADIATION HYDRODYNAMICS 



^+,V.v = 0, 
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— pV • V — ilTKpB + CKeE, 



(4) 



and 



c^Dt\p 



) 




XfF. 



(5) 



Here the convective derivative D/Dt = d/dt + v • V. The dependent quantities p, e, v, and p are 
the material mass density, energy density, velocity, and scalar isotropic pressure respectively, while 
E, F, and P arc the total frequency-integrated radiation energy density, momentum density or flux, 
and pressure tensor, respectively. The latter are the zeroth, first, and second angular moments of 
the radiation specific intensity. The specific intensity, /(x, fi, i^, i), in general varies with spatial 
position X, viewing direction fi, frequency v, and time t. Complete description of / in a three- 
dimensional space therefore requires seven scalar coordinates. The angular moments of / appearing 
in the RHD equations above may be written 



The assumption of LTE allows the source function specifying the rate of emission of radiation 
from the gas in equations (3) and (4) to be written as the Planck function B. Equations (2) 
through (5) have been integrated over frequency, leading to the fiux mean total opacity xf, and 
the Planck mean and energy mean absorption opacities, Kp and ke- In the present paper the 
opacities are assumed to be independent of frequency, so that Kp = ke and the subscripts may 
be omitted. The total opacity x is the sum of components due to absorption, k, and scattering, 
(7, all having dimensions of inverse length. In a pure hydrogen gas above 10^ K, simple forms for 
K and a are Kramers' Law and the Thomson scattering coefficient. Equations (1) to (5) correctly 
describe the flow only in the comoving frame, where material properties are isotropic and the form 
of the material-radiation interaction terms is greatly simplified. In writing the equations, effects of 
self-gravity and magnetic fields have been ignored. The treatment of these effects in ZEUS-2D is 
independent of the algorithms for RHD described here. 

The equations of RHD may be closed by the addition of constitutive relations for the Planck 
function and opacities, an equation of state specifying the gas pressure, and an assumption about 
the relationship between the angular moments of the radiation field. Since P = ^E when the 
field is isotropic, one might choose to assume this relation holds everywhere. This is referred to 
as the Eddington approximation, and is analogous to assuming isotropy of the material particle 
distribution function in deriving the gas equation of state. The Eddington approximation implies 
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that in steady-state, equation (5) becomes 



This expression gives the correct flux in optically thick regions, where x is large, and tends to 
infinity in optically thin regions where x ^ 0- However the rate at which radiation transports 
energy is finite and, from equations (6) and (7), is limited to |F| < cE. This causality problem 
with equation (9) occurs because in optically thin regions photons travel freely, their mean free paths 
may exceed characteristic lengths in the system, and the radiation field may be anisotropic. In this 
paper we implement an extension of the Eddington approximation in which causality is preserved 
by assuming a particular form for the angular dependence of the radiation field in optically-thin 
regions. This extension is described in the next section. 



3. THE FLUX-LIMITED DIFFUSION APPROXIMATION 

FLD techniques were first used in astrophysics to solve the radiation transport equation, in 
calculations of accretion onto a neutron star (Alme & Wilson 1974). Subsequently, Levermore & 
Pomraning (1981; hereafter LP) generalized the method to approximately handle transport phe- 
nomena while preserving causality in regions where spatial variations occur over distances smaller 
than a mean free path. Pomraning (1983) included relativistic terms of order v/c in a moving 
fluid. FLD methods have since been used in pure scattering media (Melia & Zylstra 1991), and for 
fully general-relativistic calculations (Anile & Romano 1992). The key ingredient in all cases is the 
assumption that the specific intensity is a slowly- varying function of space and time. In one space 
dimension, this assumption is valid in both the optically-thick diffusion limit and the optically-thin 
free-streaming limit. One hopes it holds approximately in the intermediate regime and in multi- 
dimensions, and this may be tested by comparison with results of full-transport calculations. Given 
slow spatial and time variation, analytic relations may be obtained between the angular moments 
of the radiation field. LP showed that the radiation flux may be written in the form of Fick's law 
of diffusion, as 

F = -DVE, (10) 

with a diffusion coefficient D given by 

The dimensionlcss function A = X{E) is called the flux limiter. In this framework, the radiation 
pressure tensor may be expressed in terms of the radiation energy density via 

P = fE, (12) 

where the components of the Eddington tensor f are given by 

f = ^(l-/)l + ^(3/-l)nn, (13) 
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n = V-E/|V-E| is the unit vector in the direction of the radiation energy density gradient, and 
the dimensionless scalar function / = f{E) is called the Eddington factor. The flux limiter A and 
Eddington factor / are related through implicit constraints between the moments F and P, so that 

f = X + X'R', (14) 

where R is the dimensionless quantity R = \'VE\/{xE)- 

Equations (10) through (14) close the equations of RHD, eliminating the need to solve the 
radiation momentum equation (5), and greatly simplifying the computation. Since the angular 
distribution of the radiation field is not explicitly computed, no formal solution of the transfer 
equation is needed. An important choice which must be made to achieve this simplification is the 
expression to be used for the flux limiter A. Many expressions are possible which preserve causality, 
and are consistent with the assumption of smoothness in the radiation field. These are distinguished 
by different assumptions regarding the details of the angular dependence of the specific intensity. 
In one example, LP chose to constrain the normalized specific intensity ip = I/{cE) by 

1^ + J2.VV; = 0, (15) 
c ot 

resulting in the fiux limiter 

Physically, the derivation of this result follows an analysis similar to the Chapman-Enskog approach 
for developing solutions to the Boltzmann equation in the kinetic theory of gases. In a second 
example of a choice of constraint on the anisotropy of the radiation field, Minerbo (1978) assumed 
a piecewise linear variation of the specific intensity with angle, and found the fiux limiter 



urn- 2/(3 + V9 + 12/^2) ff0<i?<3/2 

{1 + R + VT+2R)-' if3/2<i?<oo ^ ^ 

In the optically thin limit i? — ^ oo, both the LP and Minerbo flux limiters give to first order in R~^ 

hm X{R) = i, (18) 

it— ►OO It 

SO the magnitude of the fiux approaches |F| = c|V-B|/(x-R) = cE, which obeys the causality 
constraint. In the optically thick or diffusion limit — ^ 0, both examples give to first order in R 

limA(iJ) = i, (19) 

SO the fiux takes the value given by equation (9). For intermediate values of R, these two forms 
for the flux limiter differ substantially, as shown in figure 1. Thus a major uncertainty in FLD 
calculations is the appropriate choice for the form of the flux limiter. Some alternatives are described 
by Levermore (1984). The LP form has been used in studies of accretion onto black holes (Eggum, 
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Coroniti, & Katz 1988) and protostellar collapse (Bodenheimer et al. 1990). Dynamics of radiation- 
dominated accretion disks computed using the LP and Minerbo forms were compared by Turner 
& Stone (2001). Kley (1989) derived a flux limiter appropriate for accretion disk boundary layers, 
using exact solutions of the transfer equation for spherically-symmetric stellar atmospheres. This 
limiter is quite similar to the Minerbo form over the whole range of R. A flux limiter appropriate 
for a scattering medium was derived and used in accretion disk corona models by Melia Sz Zylstra 
(1991). The LP and Minerbo flux limiters have been widely used, and bracket the range of values of 
most other flux limiters. Both are available in the ZEUS-2D implementation of FLD described here. 
Since the appropriate form of the flux limiter is not always clear beforehand, in some applications 
new limiters may have to be derived in order to achieve accurate results. Comparisons between 
FLD and full transport results remain important. 

In summary, applying the FLD approximation involves computing the flux limiter from the 
value of R in each zone, and then using equations (10) to (14) to obtain the radiation flux and 
radiation pressure from the radiation energy density. 



4. NUMERICAL METHOD 

As with the MHD equations in ZEUS-2D, the dynamical equations for the radiation are 
operator-split into source and transport terms. In the source step, the radiation and material 
energy densities are updated using finite difference approximations to 

dE 

- = -V.F (20) 

and separately, 

dE 

= - Vv : P + AttkB - ckE (21) 

ot 

and 

de 

— = -pV • V - AttkB + ckE. (22) 
at 

In the transport step, an integral formulation is used to generate a conservative differencing scheme 
for the advection terms, 

^ [ EdV = - (f E{v- Vg) • dS, (23) 

dt Jv JdV 

where Vg is an arbitrary coordinate velocity which allows grid motion. Equations (20) to (23) 
are evolved in the radiation module, as is the source term due to radiative acceleration in the gas 
momentum equation (2). The rest of the set of coupled equations is solved using the existing hy- 
drodynamic algorithms. Centering of the radiation variables on the hydrodynamic grid is described 
in §4.1, update of the radiation flux term in §4.2 and the remaining source terms in §4.3, treatment 
of the compression terms in the case of strictly scattering opacity in §4.4, and the transport step 
in §4.5. Stability of the method and the choice of timestep are discussed in §4.6. 
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Fig. 1. — Dependence of the flux-limiter (upper panel) and Eddington factor (lower panel) on 
optical thickness parameter R = \VE\/{xE). Levermore-Pomraning values (eq. [16]) are shown 
by solid curves, Minerbo (eq. [17]) by dotted curves. Eddington factors are computed from the 
fiux-limiters using equation (14). 
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4.1. Centering of Variables 

The dependent variables for the radiation are discretized onto the computational mesh as in 
ZEUS-2D (Stone &; Norman 1992a, hereafter SN). Scalars, and tensors of even rank, are zone- 
centered, while tensors of odd rank are face-centered. The mesh is labeled by coordinate vectors Xi 
and X2, with the 3-direction taken to be orthogonal to the computational plane. The mesh whose 
grid lines mark zone edges is labeled by coordinates xlai and x2aj, while the mesh whose lines 
intersect at zone centers is labeled by coordinates xlbi and x2bj, as shown in figure 2. On this 
mesh, the radiation variables are discretized according to 



E{xi, 


.X2) — 


> E{xlbi,x2bj) = 


Fi,j 


Fi{xi, 


.X2) — 


> Fl{xlai,x2bj) 


= FU,, 




.X2) — 


> F2{xlbi,x2aj) 


— F2- ■ 


-Pii(a;i. 


.X2) — 


> Pll{xlbi,x2bj] 


) = PlUj 




,X2) 


> P22{xlbi,x2bj] 


) = P22,^j 


P?.?,{xi; 


.X2) 


> P33{xlbi,x2bj] 


) = P33,j 


Pi2{xi., 


.X2) 


> P12{xlbi,x2bj] 


) = P12,j 



The components of the Eddington tensor f are centered identically to the components of P. The 
opacity, absorption and scattering coefficients, and Planck function are all zone-centered. The 
radiation diffusion coefficient D is computed on zone 1-faces (Dljj) and 2-faces (Z)2jj) along with 
the components of the flux. 



4.2. Radiation Flux Divergence Term 

Under the FLD approximation, the flux divergence term in the radiation energy equation (20) 
is a function of the gradient in radiation energy density. When implicitly differenced, it therefore 
depends on the time-advanced values of E in adjacent grid zones, and may be updated by a matrix 
inversion. To simplify the form of the matrix, we operator-split this from the remaining source 
terms, and integrate 

dE 

— = V.(DVE) (24) 

with D given by equation (11). The solution method is described here for a uniformly-spaced square 
Cartesian grid. Extension to the curvilinear coordinates used in ZEUS-2D is straightforward. The 
implicit difference equations to be solved are of the form 

D2l^^,{E-+l, - E-+') - D2l^{E-+' - E^fi,), (25) 

where Ax is the grid spacing, and the spatial indices span ranges i = 1. ..N and j = 1 . . . M. 
Radiation energy densities are to be determined at timestep n+1. The difference equations are kept 
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Fig. 2. — Centering of the radiation hydrodynamic variables in the ZEUS-2D code. Symbols are 
explained in §4.1. 
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linear by using diffusion coefficients from step n. Tfiis is a good approximation when |V£^|/(x£^) 
changes httle each timestep, but may be less appropriate when the radiation field varies rapidly with 
time in regions of intermediate or low optical depth. In section 5.5 it is shown that a propagating 
radiation front can be more accurately evolved in this situation when the diffusion coefficients are 
recalculated every diffusion substep. 

Considering E^j'^ as a single vector of length A'^ x M, the updated values may be determined 
by solving an equation with a matrix of {N x M)^ elements, using a sparse-matrix technique such 
as incomplete Cholesky - conjugate gradient (Press et al. 1992). These methods have proven useful 
for particular problems with fixed boundary conditions, but the form of the matrix and the sparse- 
matrix method required vary with the boundary conditions. Solutions of equation (25) may also 
be obtained by combining simple solvers with multi-grid acceleration techniques (Hackbusch 1985). 
However, here we choose an alternating-direction-implicit (ADI) method on a single grid. This 
may be less efficient than the best sparse-matrix and multi-grid algorithms, though the number 
of operations required is proportional to N x M. In ADI, successive approximations to £'"+^ are 
computed by advancing towards a u>-stationary state the equation 

= - -'^^^ + [V • {DVE}]^;^' . (26) 

The new variable w may be thought of as the pseudo-time. Each tu-step is split into two parts. In 
the first, the update is lu-implicit along the 1-direction, and explicit along the 2-direction. In the 
second, the differencing schemes for the two axes are exchanged. Labeling the approximate value 
of E"^^^ at pseudo-timestep m by E^^, the difference equations solved are 



-,m+h ™ Aw . Aw 



E ^ -Er, = —— (Ef, - E 2 ) + 



+ ^2i:,+i(^S+i - ^S) - - ^™_i)} (27) 

on sweeps which are implicit along the 1-direction, and 

i.i i.i n\4.\ 1,3 i.i 



1,3 1,3 2AV ''^ ''^ ' 2(Aa;)2 

Di^i,(iCj - - ^1^.(47' - Ctj) 

+ ^2i:,+i(£;i^+\ - E^^^) - D2l^{E^;^ - E^^t\)} (28) 

on sweeps which are implicit along the 2-direction. Each sweep involves solving a tridiagonal matrix 
equation. For example, for boundary condition VE = 0, the matrix equation solved on the jth 



- 11 - 



sweep implicit along the 1-direction is 

" 13 - hDlij -hDl2,j 

-hDhj (5 -hDhj 

-hDhj P -hDU,j 

— /i-Dljv-ij P —hDl]\fj 

-hDlNj (3 - hDlN+i,j 

where h = 2^^, /? = 1 + ^ + h{Dli+ij + -Dljj), and the elements of vector b are 

bi = [l- h{D2,^j+^ + D2,^j)]E^^ + /iL>2,,,+i£;-.+i + hD2,^jE^^_^ + ^El^. (30) 

Advancing a pseudo-timestep involves solving M such equations in N x N matrices for the 1- 
implicit sweeps, and N similar M x M matrix equations for the 2-implicit sweeps. A standard 
method involving LU decomposition and forward- and back-substitution is used. In the case of 
periodic boundaries, the lower left and upper right corner elements of the matrices are non-zero, 
and this condition is dealt with using the Sherman-Morrison formula (Press et al. 1992) to compute 
a correction to the solution of the tridiagonal part. The diffusion coefficients are held fixed during 
the entire real timestep. Numerical accuracy of the solution is improved by scaling the lengths, 
times, and diffusion coefficients in equation (29) by a common factor chosen so the diagonal matrix 
elements are near unity. 

To bring the solution into pseudo-time steady-state rapidly on different size scales, the length 
of the pseudo-timestep Aw is increased exponentially with iteration number m= 1. .. W, so that 

A^j • 

where the initial pscudo-timcstcp Awq is one-quarter of the square of the grid spacing, and the 
final pseudo-timestep Awi is one-quarter of the square of the total grid size (Black &: Bodenheimer 
1975). In cases where the number of zones is large, the rescaling of lengths and times applied to 
equation (29) may not be sufficient to keep the condition numbers of the matrices near unity for 
all pscudo-timesteps, and the reduced accuracy of sohitions obtained using LU decomposition may 
lead to slow convergence. Under these conditions it might be helpful to iteratively improve the 
initial solution, or to use an alternative decomposition such as QR or Householder (Press et al. 
1992). For the test problems discussed in section 5, reducing the condition numbers of the matrices 
by adjusting the rescaling factor did not speed convergence. 

The error in the time-advanced solution obtained from the series of w-steps is checked 

by substituting in the original difference equation (25). The residual is normalized against the 
ratio of the time-centered E to the hydrodynamical timestep. If the biggest error on the grid is 
too large, the method is applied again with larger W and more closely-spaced pseudo-timesteps. 



E 



En,' 



(29) 
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If a sufficiently small error is not reached after many pseudo-timesteps, then the real timestep for 
the flux divergence term is halved, and the term is applied twice per hydro dynamical timestep. If 
a sufficiently small error is not obtained with the real timestep finely divided, the calculation is 
terminated. On the other hand, if the error obtained is very small, then fewer pseudo-timesteps 
are tried initially in the following diffusion substep. 



4.3. Matter-Radiation Interaction and Compressive Heating Terms 

The absorption, emission, and compressive heating source terms are approximated by implicit 
difference equations, to obtain numerical stability for timesteps which may be long compared with 
the matter-radiation interaction timescale. The solution is advanced in these terms from step n to 
step n -|- 1 by solving simultaneously 

E-;' - = At {-(Vv : P)^+i + AttkB-;' - ckE^;'} (32) 



and 

el+' - e^, = At { -pll' (V • v),,, - AttkB^J' + ckE^^' } . (33) 

These depend on the values of e and E in zone i,j only, so the spatial coordinate subscripts are 
dropped in subsequent equations. The gas is assumed to be ideal, with ratio of specific heats 7, 
and equation of state = (7 — l)e'*+^. The Planck function = ^gbT^ is computed from 

the temperature T = (7 — l)/ie"+^/(7?./9"') using the density at the preceding time. Here ^ is 
the dimensionless mean particle mass, TZ the gas constant, and cts the Boltzmann constant. The 
rate (Vv : P)"+i at which work is done on the radiation by the flow is 



{Vvn/11 + VV22/22 + Vv33(l - /II - /22) + ((VV12) + (Vv2i))/12} E"+\ 



(34) 



The Eddington tensor components /II, /22, and /12, the fluid velocity v, and the absorption 
opacity k are evaluated at timestep n, and the components of the velocity gradient tensor are 
as in Stone, Mihalas, & Norman (1992) Appendix A. With these choices, the only unknowns in 
equations (32) and (33) are e"+^ and Eliminating E^'^^ between the two equations yields a 

quartic polynomial with e"+^ as one root. With x in place of e"+^, the quartic is 

4 (l + a4)(l + a2 + a3) (1 + 02 + 03)6" + 02^" _ „ 

01(1+03) 01(1 + 03) 



where 



M7-I) 



4 



«i = 4«^b|^^^^| ^t, (36) 

02 = cnAt, (37) 
(Vv : P)"+^ 

«3 = ^ At, (38) 
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and 

04 = (7-l)(V- v)At. (39) 

The quartic typically has a single positive real root when the timestep satisfies the Courant condition 
of §4.6. Writing the coefficient of x in equation (35) as ci and the constant coefficient as cq, this 
single root lies between zero, and the smaller of |co/ci| and jcol^/'^. The updated gas energy density 
is obtained by finding the root using Newton-Raphson iteration on this interval, with bisection 
when the Newton-Raphson method fails. The updated radiation energy density is determined by 
substituting the updated gas energy density in equation (32). 



4.4. The Case of Zero Absorption Opacity 



When the absorption opacity k is zero and the material-radiation interaction terms vanish, 
equations (32) and (33) are no longer coupled. In this instance, numerical energy conservation may 
be improved by using time-centered values for the gas and radiation pressures. In this differencing 
scheme, the time-advanced gas energy density e is 



1 - igAi 



(40) 



with g = (7 — 1)(V • v)". For the radiation energy density update, q is replaced by the quantity 
multiplying in equation (34). 



4.5. Transport Step 

Advection of the radiation energy proceeds exactly as described by SN for the hydrodynamic 
variables. Having expressed the advection terms in integral form in equation (23), we apply the 
conservative differencing scheme used for the hydrodynamic variables. In this scheme, the fiux of 
advected radiation energy is computed from the mass flux in order to reduce the relative numerical 
diffusion of radiation with respect to the gas. The flux across every zone interface is computed 
using an interpolation method which may be selected as either donor cell, van Leer, or piecewise 
parabolic. These fluxes are then used to update the radiation energy density in a directionally split 
fashion. 

Thus, the flux of radiation energy density along the 1-direction is constructed using 

J^lj = iE/d)l,^^Ml^g2a'^^'^g31a';^'^dvl2a], (41) 

where the mass fluxes in the 1-dircction, M^j, arc defined by equation [55] in SN. These fluxes are 
used to update the radiation energy under 1-transport via 

E]'+^dvlla^+^dvl2a] - E^jdvlla2dvl2a] = -At(J^^^.+i - J^lj). (42) 
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Advection fluxes of radiation energy density in the 2-direction are computed by 

^2 . = {E/d)l^,^^My2>lVldxlb^g?,2a]'''^ , (43) 

where the mass fluxes in the 2-direction, Mf-, are defined by equation [56] in SN. These fluxes are 
then used to update the radiation energy due to advection in the 2-direction via 

E^f^dvlla2dvl2a]+'^ - El'jdvllafdvl2a] = -AtiTfj+i - Tfj). (44) 

To reduce the systematic effects of the directional splitting discussed in section 5.3, the order in 
which the two advection substeps are applied is reversed each timestep. 

4.6. Stability and the Timestep Criterion 

Since the advection terms in the RHD equations are treated time-explicitly, the calculations 
may be numerically unstable unless the timestep is less than the time for waves or advection to 
carry energy across a grid zone (Richtmyer & Morton 1957). The wave family with the largest 
group speed may be either adiabatic sound waves or radiation acoustic waves, depending on the 
ratio of gas and radiation energy densities. The sound speed for purposes of computing the timestep 
is therefore chosen to be 

1/2 

(45) 

where Ptot is the sum of the gas pressure and the largest component of the radiation pressure 
tensor. The remainder of the calculation of the timestep proceeds as in the hydrodynamic portion 
of ZEUS-2D (SN). Though resulting steps may be longer than the time for radiation to diffuse 
across a zone the implicit differencing of the radiation diffusion term (§4.2) ensures stability 

of the method. 

5. TEST CALCULATIONS 

The problems used to test the radiation module fall into three categories. Calculations outlined 
in the first three sections below test in isolation the heating and cooling terms (§5.1), the radiation 
flux divergence term (§5.2), and the transport terms (§5.3). Linear RHD wave tests involving all the 
radiation terms are described in §5.4. Calculations which exercise radiation terms in the non-linear 

regime are propagation of optically-thin radiation fronts (§5.5), structure of radiating gas pressure 
dominated shocks of moderate optical depth (§5.6), and formation of a steady optically-thick shock 
in which radiation pressure dominates downstream (§5.7). 



c., = 



max(7, -) 

P 
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5.1. Heating and Cooling Terms 

In a static uniform absorbing fluid initially out of thermal balance, an analytic solution for the 
time evolution of the gas energy may be obtained for the case where radiation energy dominates 
the total. The radiation energy emitted or absorbed in reaching equilibrium is a small fraction of 
the initial value, so radiation energy density E is assumed constant. The time evolution of the gas 
energy density e is obtained by solving the ordinary differential equation 

— = CKE-AiTKB{e). (46) 

Solutions are plotted together with the numerical results in figure 3, for a density of 10^^ g cm~^, 
opacity 4 x 10~® cm~^, mean dimcnsionless mass per particle 0.6, and index in the equation of 
state 7 = 5/3. Two cases are shown. In the first, the initial thermal energy density is less than 
the equilibrium value. The agreement of the computed result with the analytic solution in this 
case is excellent. In the second case, e is initially larger than the equilibrium value, lags the 
expected decline for the first few steps, and reaches equilibrium in the correct time. These results 
are consistent with the use of a fully-implicit method and a timestep of a few times 10~^^ seconds, 
much longer than the cooling time. When the timestep is held to less than the cooling time, the 
agreement between the calculation and the analytic solution is improved, as shown by the dotted 
line in figure 3. The sum of thermal and radiation energies is conserved in these three cases to 
about the floating-point precision. 



5.2. Flux Divergence Term 

An analytic solution of the diffusion equation (24) on the unit square, with periodic boundaries 
and unit diffusion coefficient, is (Lamb 1995) 

E{x, y,t) = 2 + e"^'^'^* sin 27rx sin 27r7/. (47) 

This function at t = is discretized onto a grid with 100 zones along each side, and the evolution 
followed using the radiation diffusion term alone. The timestep chosen is a factor of a hundred 
longer than the time for radiation to diffuse across a zone. The largest zone-by-zone difference 
between the numerical and analytic solutions during the evolution towards equilibrium is 5.3%, and 
the largest root-mean-square difference is 2.5%. When the timestep is shortened by a factor ten, 
the corresponding residuals are 0.7% and 0.3%. The implicit numerical solution relaxes towards 
equilibrium more slowly than the analytic solution in each case. 



5.3. Transport Terms 



When the directionally-split advection substeps described in section 4.5 are applied in the 
same order each timestep, the splitting leads to a systematic delay in transport along the first axis. 
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Logio[ time / s ] 

Fig. 3. — Results from calcuiations of approach to tiiermal equiiibrium are sliown by crosses every 
timestep. Starting value of the thermal energy density is lO^'^ erg cm~^ for the upper set of crosses, 
and 10^ erg cm~^ for the lower set. The initial radiation energy density E is 10^^ erg cm~^ in both 
cases. Corresponding analytic solutions assuming constant E are indicated by solid curves. Results 
from a third calculation, with initial energy densities matching the upper case, and timestep fixed 
at 10~^^ seconds, are shown by a dotted line almost coincident with the upper solid curve. 
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and an advance in transport along the second axis. As a test of this effect, the source terms were 
switched off, and a circular region of increased radiation energy density initially at the center of the 
domain was carried twice diagonally across a grid of 32^ zones with periodic boundaries. About 
2000 timesteps were taken, and the advection substeps were applied always in the order 1, 2. At the 
end of this calculation, the centroid of the region lay off the grid diagonal by 0.02 zones along each 
axis. When the calculation was repeated with the order of the substeps reversed each timestep, 
the centroid remained within 3 x 10~'^ zones of the diagonal. We have chosen to reverse the order 
of advection each timestep in all further calculations discussed in this article, in order to minimize 
the effects of the directional splitting of the transport terms. 



5.4. Linear RHD Waves 

Waves in a radiating fluid may be damped by escape of photons. Damping is often rapid 
for oscillations with frequency near the rate at which material cools by radiating, and also for 
oscillations with wavelength near the distance radiation diffuses in a wave period. Mihalas & 
Mihalas (1983) obtained a linear dispersion relation for radiation and acoustic plane waves, driven 
by a boundary into a uniform medium with opacity due to absorption. They made the Eddington 
approximation, and assumed a uniform background in thermal equilibrium. The coefficients of the 
dispersion relation may be specified by two dimensionlcss parameters. The Boltzmann number 
Bo = 47Cae/(c£') is a ratio of typical rates of energy transport due to material advection and 
radiation in the background state, and r = E/{4rye) fixes the ratio of the energy densities. In these 
definitions, 7 is the index in the equation of state and Ca is the adiabatic sound speed. For the 
tests discussed below. Bo = 10~'^ and r = 0.1. Three frequencies are considered, selected for optical 
depth per wavelength t\ = 10^'^, unity, and 1000. The first of these frequencies is near the radiative 
cooling rate, and the last mode has wavelength near the radiative diffusion distance. In the test 
calculations, absorption opacity is set to 0.4 cm^ g~^ times the density. The domain is initially 
stationary and uniform, and the chosen Bo and r correspond to density 3.216 x 10"^ g cm~^, 
gas energy density 26020 erg cm~^, and radiation energy density 17340 erg cm~^. A region one 
wavelength long at the left of the domain is updated each timestep with the time and space variation 
found from the linear analysis. A wave propagates to the right from this driving region. Late-time 
results for waves of the three optical depths are shown in figures 4 to 6. There are ten zones per 
wavelength in each case. Timesteps were chosen according to the Courant condition of §4.6 for 
the calculations in the left-hand panel of each figure. The resulting steps, though having Courant 
numbers of 0.5, are longer than the time {Ax)'^/D for radiation to diffuse across a zone. The 
fully-implicit method used remains stable for these long timesteps, but loses accuracy, and the 
waves are damped more strongly than expected. Similar damping occurs when the Eddington 
approximation is made, indicating the effect is not due to the explicit differencing of the radiation 
diffusion coefficient in equation (25). In the right-hand panels of figures 4 to 6 are results from 
calculations in which the timesteps were reduced so the method adequately reproduced the linear 
analytic results. The wavelengths obtained are one to two percent shorter than the predicted 
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wavelengths. Similar agreement with the hnear predictions was observed for Bo = 10 ^ and 
r = 1000. 

5.5. Propagating Radiation Fronts in Optically Thin Media 

In optically-thin regions, changes in radiation energy density propagate at the speed of light. 
To examine how well the method follows such changes, a domain 1 cm long is placed in ther- 
mal equilibrium with a radiation energy density of 10"^^ erg cm~^. The mass density is set to 
0.025 g cm~^ and the opacity to 0.4 cm^ g~^, so optical depth across the domain is 0.01. One hun- 
dred grid zones are used. At time t = 0, the radiation energy density in zones to the left of a; = 0.1 
is raised to unity, and the increase is allowed to propagate to the right. The left-hand boundary 
condition is inflow, the right-hand outflow. Results using the LP flux limiter and scattering opacity 
are shown in figure 7. The front moves at approximately the speed of light. It is spread over a 
large number of zones, indicating the method is quite diffusive in this situation. For problems in 
which details of such fronts are important, it may be better to solve the time-dependent transfer 
equation using steps set by the radiation flow time. With timesteps longer than Ax/(2c), the front 
is wider still, however the method remains stable. The width is less when the diffusion coefficients 
are recalculated every diffusion substep, rather than once per hydrodynamic timcstep. Where the 
leading edge of the front reaches the end of the grid, a region of uniform lower radiation energy 
density is seen. The effect is greatest when the timestep is long and the front is broad, as in the 
curve shown by dots on the bottom panel of figure 7. Radiation energy accumulates in this region 
because the boundary condition VE = forces the radiative flux to zero at the boundary. Other 
conditions on radiation energy density which may be appropriate across outflow boundaries are 
fixed values of E, and V • F = 0. 

The off-diagonal components of the Eddington tensor (equation (13)) may be non-zero in 
radiation fronts inclined with respect to the coordinate axes. As a test of the method in such a 
situation, the domain used in the one-dimensional radiation front calculation is expanded to a 1 cm 
square of 100 x 100 zones. The high radiation energy density is placed initially in zones within 
0.1 cm of grid center. The boundaries are made periodic in both directions, and the timestep 
is fixed at Ax/(2c). The distribution of radiation after 10^^^ seconds is shown in figure 8. The 
lowest contours on the figure are non-circular. After the leading edge of the front reaches the 
midpoints of the grid boundaries, the gradient in radiation energy density there is directed along 
the boundaries, rather than away from grid center. Since in the FLD approximation the flux is 
assumed parallel to —VE, an incorrect result is obtained in these regions. A similar error will 
result in FLD calculations whenever inclined radiation fronts meet. After the front has passed the 
corners of the grid, at 5 x lO"-*^^ sec, total energy differs from that at the start of the calculation 
by about one part in 10^^, indicating that energy is conserved to high precision in this instance. 
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Fig. 4. — An acoustic oscillation driven in the region left of the vertical dashed line is damped 

by radiative cooling as it travels to the right. Optical depth per wavelength is 10~^. Dotted lines 
indicate the linear analytic solution, solid lines the numerical solutions. The quantities plotted 
are the relative perturbations in (top to bottom) density, thermal energy density, radiation energy 
density, and velocity. There are ten grid zones per wavelength A. For the left panel, the timestep 
is set by the Courant condition, and is longer than the time for radiation to diffuse across a zone 
by a factor 1.5 x 10^. For the right panel, the timestep was reduced by a factor ten. 
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Fig. 5. — Same as figure 4, but for a slower oscillation such that the optical depth per wavelength 

is unity, and damping by radiation is slight. There are ten grid zones per wavelength A. For the 
left panel, the timestep is set by the Courant condition, and is longer than the time for radiation 
to diffuse across a zone by a factor 1.5 x 10^. For the right panel, the timestep was reduced by a 
factor 1000. 
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Fig. 6. — Same as figure 4, but for a still slower oscillation such that the optical depth per wavelength 
is 1000, and damping occurs by radiative diffusion. There are ten grid zones per wavelength A. 
For the left panel, the timestep is set by the Courant condition, and is longer than the time for 
radiation to diffuse across a zone by a factor 15. For the right panel, the timestep was reduced by 
a factor 100. 
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Fig. 7. — Radiation front propagating across a box of optical depth 0.01 is shown at three different 
times. The solution obtained with timcstcp A.t/(2c) is shown by dashed curves. Results with 
timestep a factor ten longer are indicated by dots, and those with timestep a factor ten shorter by 
solid curves. Vertical lines indicate the expected positions of the front. 
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Fig. 8. — Radiation front expanding from near the center of a square of optical depth 0.01 per 
side. Region where radiation energy is initially enhanced is marked by a heavy contour. Expected 
location of the front after 10~^^ sec is shown by dashed circle. Calculated distribution of radiation 
energy density at this time is shown by light contours at 1%, 25%, 50%, and 75% of the maximum. 



-24- 



5.6. Subcritical and Supercritical Shocks 

The structure of a shock in radiating fluid is sometimes altered by photons propagating up- 
stream to preheat approaching material. If material reaching the shock front remains cooler than 
the shocked gas, the shock is termed subcritical, and the main effect is a slightly raised post-shock 
temperature, which declines downstream as radiation is emitted upstream. In stronger, supercritical 
shocks, upstream material is heated almost to the temperature of the post-shock gas, and material 
immediately upstream is in thermal equilibrium with the radiation emitted through the front. Gas 
reaching the front then overshoots the flnal temperature significantly, and cools to the post-shock 
temperature by radiating from a layer with optical depth less than unity. In this section we compare 
structures calculated for subcritical and supercritical shocks against analytic estimates by Zel'dovich 
& Raizer (1967), and numerical results obtained by Sincell, Gehmeyr, & Mihalas (1999) using a 
radiation hydrodynamics code with adaptive refinement of a one-dimensional grid. The problem is 
set up as by Sincell et al. (1999). The gas is given absorption opacity 0.4 cm^ g^^, mean molecular 
weight 0.5, and uniform density 7.78 x 10"^*^ g cm~^. Gas temperature 10 + 75a;/(7 x lO^'' cm) K 
increases linearly with distance x from the origin, and the radiation energy density is chosen for 
equilibrium. Grid spacing is set to 1.4 x 10^ cm, so initially the optical depth per zone is 0.044. A 
piston moves into the gas from the origin at uniform speed, and shocked material accumulates in 
front of the piston. 

Results for a subcritical shock, with piston speed 6 km s~^, are shown in figure 9. The time is 
10^ seconds. Temperature declines from the post-shock maximum more rapidly than expected, and 
the fiux is too large in several zones downstream from the front. These effects may be due to an 
incorrect angular distribution of specific intensity. They are similar when the numerical resolution is 
doubled, suggesting they are not due to the finite optical depth of the front imposed by the artificial 
hydrodynamic viscosity. Sincell et al. (1999) found, by directly solving the transfer equation, that 
the Eddington factor near the shock front is slightly less than | because the emitting layer has 
greatest line-of-sight optical depth when seen at grazing angles. In the FLD approximation, angular 
variation is assumed to be either smooth (LP, eq. [15]) or stepwise (Minerbo), and this feature is 
missed. Upstream, flux is higher with the LP than with the Minerbo flux limiter. As may be 
seen from figure 1, the flux is reduced below the Eddington value at larger optical depths with the 
Minerbo limiter. 

The structure calculated for a supercritical shock with piston speed 16 km s""^ is shown in 
figure 10. The time is again 10^ seconds. The temperature spike at r = is about one-fifth as 
tall as found by Sincell et al. (1999), and about five times as wide, due to the limited spatial 
resolution. The peak temperature in the spike is 314 K above the downstream temperature Ti, 
whereas Zel'dovich &: Raizer estimate that the spike amplitude is — l)7i = 1619 K. With 
spatial resolution doubled, the spike is narrower and its height is larger at 506 K. The feature has 
approximately the expected degree of asymmetry, extending to a larger optical depth downstream 
than upstream. The flux in the equilibrium preheated region is slightly overestimated, but away 
from the front varies with optical depth in a manner close to that predicted by Zel'dovich & Raizer 
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Fig. 9. — Temperature (top), flux (middle), and Eddington factor (bottom) in a subcritical shock. 
Optical depth r is measured from the shock front. Parameters are the same as for Sincell et al. 
(1999) figures 4 and 8. Results obtained using the LP flux limiter are shown by squares, Mincrbo 
limiter by crosses. Approximate analytic solutions from Zel'dovich &: Raizer (1967) are indicated 
by solid lines. Flux is normalized to the value at the front in the analytic solutions. Solid line in 
the bottom panel indicates Eddington factor i. 
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(1967). The flux upstream from the equihbrium region depends on the choice of hmiter, with the 
Minerbo Hmiter producing smaUer flux and lower Eddington factor, as in the subcritical case. 

In both subcritical and supercritical shock calculations, the outcome is independent of the 
timestep for steps smaller than about ten zone radiation diffusion times . We conclude from 
these results that on a uniformly spaced grid, the FLD approximation is a poor choice for deter- 
mining detailed structure of shocks of intermediate optical depth. It may be sufficiently accurate 
for multi-dimensional calculations in which such shocks must be followed approximately as part of 
a more complex flow. 



5.7. Radiation-Dominated Shock 

If material upstream and downstream from a radiating shock have sufficient optical depth, 
radiation is trapped near the interface, diffusing a limited distance upstream. In this section we 
examine whether the code adequately matches the jump conditions and thickness of such a shock. 
Parameters are chosen so the energy density is mostly in the gas upstream and in the radiation 
downstream. Opacity is assumed to be due to absorption as in §5.6, and material and radiation 
are assumed to be in thermal equilibrium far from the shock. Initially, the left of the grid is set to 
density 0.01 g cm~^, temperature 10^ K, and speed 10^ cm s~^, yielding a Mach number of 658. The 
right is set to density 0.0685847 g cm~^, temperature 4.239x 10^ K, and speed 1.458 x 10^ cm s~^, as 
computed from the jump conditions. The boundary conditions are inflow on the left, and outflow on 
the right. After a brief transient, a steady shock is established. The situation after 50 flow crossing 
times is shown in flgure 11. The approximate expected thickness of the shock is the distance / for 
which the time to diffuse upstream, /D, matches the time to sweep downstream, l/ui (Mihalas 
k, Mihalas 1984). Here D is the radiation diffusion coefficient and ui the upstream gas speed. The 
shock thickness obtained in the calculation is consistent with the resulting estimate of 7500 cm. 
The values downstream in the calculation differ by a maximum of 1.0% from those predicted using 
the jump conditions. When zone spacing is decreased a factor four, the discrepancies are reduced to 
0.2%, indicating that mass, momentum, and energy are adequately conserved in material passing 
through the shock. 



6. SUMMARY 

We have described a set of subroutines which may be used with the ZEUS-2D magneto- 
hydrodynamics code for radiation MHD calculations in the flux-limited diffusion approximation. 

Results were obtained typically in three to ten times more floating-point operations per timestep 
than required for similar calculations without radiation, with the greater number needed in the 
optically-thin tests. The primary limitation of flux-limited diffusion is that the flux is not evolved 
as an independent variable, but is assumed anti-parallel to the gradient in radiation energy density. 
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Fig. 10. — Gas and radiation temperatures (a), closeup of the post-shock temperature spike (b), 
flux (c), and Eddington factor (d) versus optical depth in a supercritical shock, with LP flux 
limiter. These may be compared with Sinccll ct al. (1999) figures 7 and 9. Radiation temperature 
in (a) is shown by a sohd line, gas temperature by points. An approximate analytic solution from 
Zel'dovich &: Raizcr (1967) is indicated by a dashed curve, and horizontal dashed line indicates the 
predicted height for the temperature spike. Radiation temperatures in (b) are shown by squares, 
gas temperatures by triangles. Flux in (c) is normalized to the value at the front in the analytic 
solution. Dashed line in (d) indicates Eddington factor ^. 
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Fig. 11. — Shock of Mach number 658 in an optically-thick absorbing medium. Horizontal dashed 
lines show upstream and downstream states expected from the jump conditions. Vertical dashed 
lines indicate the approximate expected extent of the preheated layer ahead of the shock. Optical 
depth between the vertical lines is about 80. Each cross marks one of the 100 grid zones. Quantities 
plotted, from top to bottom, are the density, thermal and radiation energy densities, and velocity, 
in cgs units. 
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This is incorrect in flows where shadows are cast or incUned radiation fronts cross. The incorpora- 
tion in the flux hmiter of a specified angular dependence for the radiation field has greatest effect 
in regions of intermediate to low optical depth, and may contribute to the broadening of radiation 
fronts. Flux limiters derived for angular dependences occurring in particular applications may be 
useful in obtaining more accurate solutions. 

The limitations of FLD revealed in the tests mean that care must be taken in applying the 
module to problems with highly anisotropic radiation fields in optically thin regions. In such 
cases, full transport methods may be more appropriate. However, in optically thick flows, FLD 
is more eflficient. The module has proven useful in studies of the internal dynamics of optically 
thick, radiation-dominated accretion disks (Agol et al. 2001; Turner &; Stone 2001). Although the 
method remains stable for long timesteps, it shares with other implicit schemes the property that 
fluctuations with periods less than the timestep may be rapidly and unphysically damped, as seen 
in the linear wave and sub- and supercritical shock tests. When the diffusion timestep constraint 
of §4.6 is violated, it is advisable to check that results are independent of the timestep. 

The radiation terms included in the equations solved are those of order unity in vjc in at 
least one of the free-streaming, static diffusion, and dynamic diffusion regimes (Mihalas h Mihalas 
1984; Stone, Mihalas, h Norman 1992). Adding higher-order terms would allow study of special 
relativistic effects such as radiation drag and photon viscosity. In the implementation described 
here, the radiation source terms are operator-split, with the flux divergence term evolved separately 
from the remainder. Linear waves and shocks to which both groups of terms contribute were evolved 
adequately, but situations may be encountered where the operator splitting leads to inaccuracies. 
We have assumed the gas is in local thermodynamic equilibrium, and relaxing this assumption may 
in the future allow a great variety of additional problems to be addressed. 
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